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Abstract 


The present work extends the recently 
reported implicit analogue of MacCormack's 
earlier widely used explicit method to external 
axisymmetrlc laminar flows with strong entropy 
gradients. The details of the "numerics" of the 
implicit part are provided in a body-oriented 
coordinate system with a moving outer (shock) 
boundary during the transient part of the solu- 
tions. The limiting values of the Courant number 
are obtained when the shock boundary is treated 
explicitly. The solution algorithm outlined 
includes the treatment of the source term 
associated with the equations in weak conserva- 
tion form. From the results obtained for two 
sample problems, it becomes clear that accuracy 
of predictions is, Indeed, very good at higher 
values of the Courant number. There is a 
significant saving in overall computing time, 
depending on the Courant number used and the flow 
Reynolds number. These properties combined with 
the simplicity of programing the implicit 
analog may appeal to researchers for using 
it in the analysis of 3-D flow problems. 

Nomenclature 

constants with values less than or 
equal to unity 
skin-friction coefficient, 

(2 (9u/3n)^ 

heat -transfer coefficient, 

(2 p^/PrRe) (3h/3n)^, 

Courant number 
speed of sound, y yp/p^ 

constant in Sutherland *s law of 
viscosity, K 

nondimens ional total enthalpy, 

* . *o 

H /U„2 

nondimens ional specific enthalpy, 

h*/U*2 

finite-difference point in s-directlon 
finite-difference point in n-dlrection 
count of time steps 
freest ream Mach number 
molecular weight of mixture 
coordinate direction normal to the 

it a 

body, n /Rjj 
Prandtl number 

nondimens ional pressure, P*/p^U*^ 
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freestream Reynolds number, 

body radius normal to the body axis, 
r*/RN 

coordinate measured along the body, 

S /Rjj 

nondimens Ional temperature, T*/T* 
freestream temperature, K 
nondimens ional time, t*U*/RjJ 

freestsream velocity, m/s 
nondimensional tangential velocity, 

it it 
U /Ucc 

nondimensional normal velocity, v*/U* 
= ( u ^ + v 2)/2 
shock angle 

= Y-1 

= r+ncos0 

mesh refinement parameter 

transformed coordinates along the body 
surface 

transformed coordinates normal to the 
body 

ratio or specific heats 

shock standoff distance, 6*/RjJ 

local curvature, k /r^, also time 
step counter 
= 1+ nK 

nondimensional viscosity 

freestream viscosity, Ns/m^ 

nondimensional density, p*/p* 

freestream density, kg/m^ 

transformed time variable 
body angle 


Superscript 


* dimensional quantity 

Subscripts 


conditions at the axis of symmetry 
conditions at the wall 
conditions in the freestream 


Introduction 


past in developing computationally efficient 
methods for solving the equations of compressible 

these methods are 

the implicit time-dependent finite-difference 






1 



techniques^ which are not subject to the 
conventional stability condition of explicit 
methods.'^ However, the application of these 
techniques is frequently limited by the large 
computer time per step, their programing 

complexity, as well as severe accuracy criteria* 
These limitations increase in severity in three- 
dimensional flow analysis. In 1981, MacCormack 
presented^ an implicit analog of his earlier 
widely used explicit method.^ One of the basic 
features of this implicit analog is that it 
involves the inversion of only upper or lower 
block bidiagonal matrices as opposed to the more 
costly Inversion of block tridiagonal matrices 
needed in the existing implicit methods . ^ 

The other major advantage with the method of 
Ref. 5 is that with more complex problems, only 
the explicit part of the code Increases in 
complexity. The implicit step, which is simply 
the numerics to obtain enhanced stability with 
larger values of the Courant numbers, is not 
affected. Since it is easier to program the 
explicit part even with the increased 
complexities, the potential for this method is 
greatly enhanced, especially for the 3-D flow 
problems. Moreover, running a program fully 
explicitly can provide solutions for comparison 
and check when the implicit analog is used. 

This is an important feature for the accuracy 
check not available with the other Implicit 
methods. The implicit part in the MacCormack *8 
new method^ is merely an "add-on" to the explicit 
part.^ 

Recently, Refs. 6, 7, and 8 presented 
solutions for internal flow problems whereas 
Ref. 9 provided results for external transonic 
flows with an integral formulation by using the 
new implicit method. References 6 and 8 have 
basically used the MacCormack »s method in 
Cartesian coordinates as presented in Ref. 5. 

The results of Ref . 7 were obtained in more 
general coordinates ^(x,y), n(x,y) for a fixed 
outer boundary and Ref. 9 employed the Cartesian 
velocity vectors in the solution vector. Except 
for Ref. 8, the equations solved in these analy- 
ses were of strong conservation form. Reference 
10 has outlined a procedure for the governing 
equations which appear in "weak" conservation 
form. In this form the source terms, which are 
introduced into the equations by coordinate 
transformation and/or by turbulence modeling, 
appear outside the derivatives of the conserved 
variables. 


Analysis 

Flow Governing Equations 

The time -dependent viscous -shock-layer 
equations employed in the present analysis can be 
obtained from the unsteady Navier-Stokes 
equations by keeping terms up to second order in 
the Inverse square root of the Reynolds number in 
both the viscous and inviscid regions of the 
shock layer. These equations, when represented 
in the body-oriented coordinate system (see Fig. 
1) for a perfect gas flow at zero angle of 
attack, are expressed as^^ 


3t as 


‘IS-* 


O) 


where the vectors U, M, N, and Q may be obtained 
by dropping the species-continuity equation and 
taking the mass fraction of species and radiative 
heat flux as zero in Eq. (1) of Ref. 11, These 
vectors are also given in Ref. 12* 

The following limiting form of the governing 
equations is obtained at the axis of symmetry by 
differentiating Eq. (1) with respect to s and 
taking a limit as s ->■ 0: 


au 

o 

at 


+ 




0 


( 2 ) 


where the vectors Uq, Mq, Nq, and may 

once again be obtained from Eq. (2) of Ref. 11 or 

from Ref. 12. 


The equation of state is given by 

* * it it r. 

p » [R T /M llo^lpT (3) 

and the laminar viscosity is obtained from the 
Sutherland's law 


M = [(1+ c /T )/(T+ c /T (4) 

Transformation to Computational Plane 


The first of the two independent transforma- 
tions employed maps the physical domain into a 
rectangular region in which both the shock and 
the body are made boundary mesh lines of the 
computational region. This transformation is 


In the present analysis, MacCormack *s 
implicit analog has been extended to external 
axisymmetric laminar flows with strong entropy 
gradients. The matrices involved in the 
numerics" of the implicit part have been 
obtained in a body-oriented coordinate system 
with a moving outer boundary. The limiting 
values of the Courant number are provided when 
the shock boundary is treated explicitly. The 
method switches automatically from implicit to 
fully-explicit mode whenever the time step. At, 
satisfies the explicit stability condition. In 
general, the method becomes implicit only in 
regions where the gradients of the flow variables 
are large and a refined mesh is needed for higher 
accuracy. 


T = t, ^ E s, and n « 1 - n/6(s,t) (5) 

The transformed forms of Eqs. (1) and (2) 

are 


31 

3t 


3n 


(6) 


au 


9t 


3M 


3^ 


3N 

3n 


where 


U = U/J, M 


Me^/j 


( 7 ) 
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For n odd 


S - (U + M Hg + N n„)/J. Q - Q/J 

“o - 

% - (U/, + + N^;^„)/J. - QJJ 

Here, J is the transformation Jacobian given 
as 


, „ a(T ,C,n) ^ j, 

3(t,s,n) " ^ 


( 8 ) 


and Cs denotes 3C/3s, and so forth. 

The computation region is next mapped to 
another plane to allow higher resolution through 
the viscous layer near the surface. The second 
transformation is 



With this transformation, the final forms of 
Eqs. (6) and (7) become: 


3U ^ 3N ^ 


(10) 

n » 



3U 3M 

3N 


— ° + — 2. 
3t 3C 


(11) 


where 


U = U/J, M = M ^-/J 

R = N n-/J, Q = Q /j 

n 

and 

K = K = 

% “ ^-/j. \ = ^/j 


Here, once again, J is the transformation 
Jacobian given as 


j- / 1 \ 

9(t,C,n) \ / 


(12) 


Elements of the Numerical Integritlon Method 


Since the numerical Integration method 
has been adopted from Ref. 5, the development of 
the method presented here for a body-fitted 
coordinate system with the moving outer boundary 
will not contain the details provided there. The 
method outlined here will provide details more 
specific to the problem being analyzed. Equation 
(10) may be integrated in time by the following 
Implicit predictor-corrector set of finite- 
difference equations: 


Implicit part: 




I Implicit part : 
A 


P: < 


I - At 


* s 







AU, 


i.j 


ufy = t + 




Explicit part: 


AU^= -Ax^- 
Implicit part: 


^A M 


,k+l 


A \ 

-"l.J , ^-*^1.3 

AC An / 


(13) 


At Q. 
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And for n even 


Explicit part: 


Au; 




-At 


"k 
AC 




«k 


Implicit part: 



Tfixplicit part: 


i.j \ A? An / 


- AtQ 


■TE+r 


i.j 


Implicit part : 







where the source terms ^and ^ represent 
parts of the total source identified with the 
coordinates ^ and n, respectively, as discussed 
in Ref, 10. Based on the definition of Q given 

in Ref. 12, ^ and ^ are provided in 
Appendix A. 


It may be noticed from Eqs. (13) and (14) 
that the order of differencing has been reversed 
in the streamwise direction (while retaining the 
same order in the direction normal to the 
surface) for both the predictor and corrector 
steps between consecutive time steps. Further » 
in these equations, (A 4 /AC) and (A^yAn) are 
one-sided forward differences and (A_/AC) and 
(A_/An) are one-sided backward 

differences. ISJ and jli are matrices with 
positive eigenvalues aftd are related to the 

Jacobians A - Ofi/au) and B - ON/aS), respec- 
tively, and I is the unit matrix. The dots in 
the implicit steps indicate that the 
difference operators (A 4 . / Ac|), (A+ /IatiI), 
etc., apply also to all factors to the right. 

Both predictor and corrector consist of two 
parts : the explicit part solves the governing 

equations subject to restrictive explicit 
stability conditions ^ the implicit part removes 
these conditions by numerically transforming the 
equations into an implicit form. 

The Jacobians A and B are related to the 
Jacobians A/ A = 3M/3U and B « 3N/3U by 


® ^ ''t + r % + ® % 


(15) 


Where Ct = = 3^/38, “ 3^/3n, 

r\t * 3n/3t, Tig « 3n/3s, and tin “ 3n/3n. 


The Jacobian matrices 1 and B are provided 
in Appendix B. In the definition of these 

matrices u and v are the contravarient velocities 
obtained from 




Now, the integration method contained in the 
finite-difference Eqs. (13) and (14) can be 

simplified if the matrices A and B are 
diagonalized. Knowing the eigenvalue matrices 

Aa and Ag for A and B, respectively, these 
matrices are diagonalized as: 


1 






I 


§-1 

n 




(16) 


where and \ are the eigenvector matrices 
of %. and B and ^and Sj^ ^are the inverses of 


4 


and respectively. The complete 
definition of these matrices is given in 
Appendix C, 

The diagonal matrices Aa and Ag in Eq, 

(16), formed from the eigenvalues of X and 8, 
respectively, may be written down as: 


"A1 


0 

0 

0 


"A2 


^B1 


0 

0 

0 

a 

0 

0 


A3 


B2 


B3 


0 

0 

0 

"^A4 

0 

0 

0 


(17) 


(18) 


where 


^A1 “ ®A2 " ®A3 ° °A4 “ 


°B1 " °B2 “ °B3 “ ®B4 ° '^"8 


with 


-f 


+ (e_)‘ 


f 


g = c|/(ti A)2 + (r, )2 


The matrices X and 8 appearing in Eqs. 
(13) and (14) may now be finned by replacing the 
matrices Aa and Ag by positively-valued 
diagonal matrices Da and Dg. The matrices 

|a| and jSj are thus defined as 


“a \ 


RJ fn — 1 a 

B - S ^ S 
n B q 


(19) 


where 


- max (jA^j + I, 0.0) 
Dg - max (|Agj + Ag I, 0.0) 


( 20 ) 


and 


A 

A pA? 


1 ^ 
B pAn 




- c 


. ^ 

'A At 

An 
B At 


(21) 


time step At permitted in the n -direction, for 
example, through the Courant -Friedrichs -Lewy 
(CFL) criterion 


At < CN 


An 


V + g + 


2ym 


pRePrAn 




( 22 ) 


Here CN is less than or equal to unity for an 
explicit stable solution. 

Following Ref, 10, a procedure can also be 

outlined for obtaining the Jacobians (3^/3U) 

and (9Qg/3U) of the source term Q(“Qa + %) 

with respect to U. Justification of this 
procedure is provided in Ref. 10. 


Let = At 


and « At 


9U 

w 

3® 


(23) 


We now define the scalar matrices and 

as follows * 


|*B = '^B ^ 


Here ij)^ and (j)g are obtained from® 

(j)^ > max (q^ - Ax p^, 0.0) 
o 

't’R >. (5 r - At p„, 0.0) 


(24) 


(25) 


with 


-f 

O ' 

={" 


At , 

max — (a^ + 

max 


At . « ^ 

An ^B^ 

max 


A^) > 0.0 1 

. o.o| 


(26) 


and 


= I^ajI 

Pg - max |Igj I I 


J • 1, 2, 3, 4 


(27) 


where A^j Is the jth eigenvalue of (3^/30) 

(26) are 

obtalned“f?om 


with V ■ max|4 ^ ^ -^1 

(3 Re ’ Re Pr) 

The constants and Cg are related to 
the Courant number used for explicit stability. 
The Courant number, CN, is related to the 


max a 


max j 


Aj 


Og - max a 
max j ' J 


j - 1, 2, 3, 4 . (28) 
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and OAj defined following 

Eq. (18), whereas given by 

Eq. (21). 

The Constants and Cg appearing in 
Eq. (21) are assigned a value of 1/2 in Kef. 5. 
However, if the time step At in the expressions 

« M 

for Xa and Xg is a local minimum, a value of 
up to unity may be used for Ca and Cg to 
speed up the calculations. 

Solution Algorithm and Boundary Conditions 

As shown in expression (20) for the diagonal 
matrices Da and Dg, their elements are 
non-negative. Whenever the elements become 
negative, they are replaced by zero. This 
implies that the CFL condition, Eq. (22), is 

satisfied and matrices 111 and 111 become null 
matrices. For this case the implicit portion of 
the scheme contained in Eqs. (13) and (14) will 
be bypassed. For the flow regions where the CFL 
condition is not satisfied, the implicit parts in 
Eqs. (13) require the solution of upper-block 
bidiagonal system of equations for the predictor 
step and the solution of lower-block bidiagonal 
system of equations for the corrector step, etc. 
The integration scheme , for the case where the 
CFL criterion is not satisfied, can be illustrat- 
ed by solving the predictor part of Eq. (13) as 
detailed in the subsequent paragraphs. 

In this algorithm we replace the matrices 
$A Og (defined in Eq. (23)) by the scalar 
matrices I$a |0b| provided by Eqs. (24) 

and (25).’ ' ' ' 



then, the predictor step becomes 



an upper -bidiagonal equation and the solution can 
be obtained for each j by sweeping in the 
decreasing i direction. 

After obtaining for all i,j , then 



+ 



k 

l*j+l 


6U 




i,j+l 


(30) 


This equation is also upper-bidiagonal and is 
solved for each i by sweeping in the decreasing j 


direction. This gives for all i,j. Then 


«k *k+i 


For understanding the present method, let us 
examine the procedure for solving the block- 
bidiagonal Eq. (30) for the n-coordinate or 
j -direction. If we define 


"1 


‘"i.j * 


I An 


1 


-R.-TX 


i,j+l 


(31) 


the Eq. (30), after some matrix multiplication, 
may be written as 

where we have substituted for IbI from Eq. (19) 
and employed the relation Ugl • ^gl for the 
scalar matrix 


The integration procedure in the implicit 

part begins with the vectors ^ x, given for 
all i « 1 , 2, ..., I and j » 1, 2, .. ., J; 

•SUi,j given for all i - 2, 3, ..., I-l and j - 

2, 3, ..., J-1 and (®li,J “11 

1 = 2, 3, I-l. The quantity [B|Xj 

represents the flux of change that 'crosses tfie 
top mesh boundary. If this boundary is located 
in the far flowfleld or if the mesh is stretched 
so that At satisfies the local explicit stability 
condition (22) at the mesh points near the 
boundary, as in the case of the test problems to 
be discussed later , this flux is set equal to 
zero. Otherwise, it should be suitably specified 
from the boundary conditions. 


Following Ref. 5, the solution algorithm may 
now be summarized for each i and for j ■ J-1, 

J'"^, ..., 3, 2 in the following seven steps: 


1 ) 


1 ( An 


( 

|An 





2) X. 


(s Y w, 

I Vi.j ^ 

3) Dg - max | Ug I + Xg I, 0.0 } 


6 



4) Y, 


(I +«B 


l.j 


)I + 


1-1 


Iah 


*'b 

i.j i.j 


6) Z - Y, 


conditions are used in the ^-direction, implying 

the end flux terms jXI^ . IXI^ j 

6Uj 4 etc* to be zeto! *ln fact, in'tfiis case 

the implicit part in the ? “direction is bypassed* 

The explicit boundary conditions employed 
are no slip at surface, no surface mass transfer, 
a specified wall temperature, and pressure at the 
wall is assumed to be equal to the pressure at 
the adjacent grid point in the normal direction. 



The matrix inversion of step 4 is trivial 
because the matrix Dg is diagonal. The 
solution at grid point (i,j) is obtained at step 
5, and the flux to be used at grid point (i,j-l) 
is obtained at step 7. In computational plane 


( An )i,j+i - ( An )i,j “ An and, If global 


minimum time step 


8 used, (At)i 4+1 

IT 


(AT>i^j “At. For this case in step 1 

is obtained from 


6U 




k 

i,J+l 


6U 




i>j+l 


However, if a local minimum time step is employed 
in the computational plane (with )i,j+l ~ 

( I An I ) > Wj may be obtained from * 


(At) 


W. 

J 


6U. 


i>j (At) 


i.j 


i,j+l 




i.j+1 


For the boundary condition required at the 
solid wall boundary, the computed end flux terms 

, are saved for use as boundary condi- 
the corrector step that sweeps away from 
this boundary in the increasing j -direction. 

Using the reflection principle for the wall 
placed between the first and second grid point , 
the starting flux for the corrector step is 
obtained from 



^ _ eIbI 

1,1 ' 5 ^ 1,1 - 


1.2 ^“ l ,2 


where 


1 

0 

0 

0 


0 0 

1 0 

0 -1 

0 0 


0 

0 

0 

1 


This condition ensures that the net mass, 
tangential momentum, and energy fluxes trans- 
mitted across the wall vanish and that the net 
transverse momentum at the solid wall remains 

zero between k+1 and k+1 time steps. 


The Rankine-Hugonlot relations are employed 
for the explicit outer-boundary conditions to 
obtain flow properties immediately behind the 
shock. These relations in the body-oriented 
coordinate system are provided in Ref. 13. The 
flow conditions along the supersonic downstream 
boundary are obtained by extrapolation from the 
upstream grid points. 

Artificial Damping 

A fourth-order damping^** is used in the 
explicit part of the corrector step in Eqs. (13) 
and (14) for obtaining stable solutions over a 
large number of time steps* The following 
damping term is used in both the predictor and 
corrector steps with the implicit parts: 



This term is evaluated during step 3, of the 
solution algorithm given in the previous section, 
using the first element of the vector X. 
Accordingly, Dg in step 3 is obtained from 

Dg = max I IAj I + I + I, 0.o| 

As steady state is approached Ixj ^ I approaches 
zero and the added term vanishes • ' 

Discussion of Results 


The numerical method presented here has been 
applied to two test problems , both involving the 
analysis of viscous -shock-layer equations in the 
body-oriented coordinate system. These two 
examples are taken from Refs. 12 and 15 and 
provide fairly severe viscous-shock-layer flow 
fields for testing the present method. The main 
difference between the two test problems is that 
the first one (taken from Ref. 12) is character- 
ized by a reference Reynolds number (based on 
nose radius and freestream conditions) with a 
value of about 1.57 x 10^, whereas the second 
problem^ ^ has a reference Reynolds number around 
1.23 X 10^. The flow conditions of Ref. 15 are 
considered typical of the Jovian entry 
conditions. The reference Reynolds number 
mentioned here is related^ ^ to the mesh Reynolds 
number and provides a criterion by which the mesh 
near solid-surface boundaries may be refined. 


In the present work the Courant number 
employed in the ^-direction was always less than 
unity due to the large mesh size employed in that 
direction. Accordingly, purely-expliclt boundary 


Problem I - (Itef. 12) 

Probe geometry: 45® half -angle 

spherically-blunted cone with a nose radius 
(Rn ) of 0.222 m. 
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Jovian atmosphere: Hydrogen-helium mixture 

(0.90 H 2 + 0.10 He) under perfect-gas conditions. 


specification of boundary values for the implicit 
part near the shocks is limited by the relation 


Other flow field parameters: - 43.84, 

T«* - 145K, P«* - 1.27 X 10"“ kg/m^, 

- 4000K, Y - 1.224, Kg* - 3593.6 
J kg-* K"*, Re - 1.567 x 10^, Pr - 0.72. 

Problem II - (Ref. 15) 


Probe geometry: 44.25® half -angle sphere- 

cone with a nose radius (R^*) ot 0.352 m. 

Jovian atmosphere: Orton nominal atmosphere 

of hydrogen-helium mixture (0.895 H 2 + 0.105 
He) under perfect-gas assumption. 


Other flowfield parameters: * 43.76, 

T«* - 151.2 K, p„* - 4.966 x 10-“ kg/m^ , 

T„* » 4022,80 K, Y “ 1.217, Rg* - 3737.45 
J kg-* K-*, Re - 1.227 x 10^ , Pr - 0.72. 


Through the transformations of Eqs. (5) and 
(9), the physical domain, shown in Fig. 1, is 
transformed into a computational domain with 
equally spaced grids in both the directions, 
along and normal to the body surface. Parameter 
0 in Eq. (9) controls the amount of grid refine- 
ment desired near the wall in the physical 
domain, with values near 1 giving the largest 
amount of refinement. However, the mesh refine- 
ment is done only to the point at which the mesh 
Reynolds number reaches the order of 
unity^^>^^. At this point, diffusion and 
convection processes are equally resolved. To go 
beyond this point to smaller mesh sizes and hence 
to lower mesh Reynolds numberst , one would arrive 


at a mesh scale at which diffusion dominates. 

For such problems, more complex and time- 
consuming methods^ involving tridiagonal- 
inversion procedures should be used. In the 
present calculations, the value of 3 was chosen 
to obtain the mesh Reynolds number of order 
unity and pAx/p(An)^ kept less than 1/2 to 

avoldany possible steady-state solution 
dependence on At • This was done by reducing the 
time step near the end of the calculation. The 
damping coefficient^** e used with the explicit 
part was also reduced with the reduction in At so 
as to permit comparisons between the solutions 
having similar amounts of artificial damping. 

The damping term associated with the Implicit 
part goes to zero as steady state is approached. 


The various results presented here have been 
obtained for the values of the Courant number 
(CN) ranging from 1 to 15. The maximum value of 
the Courant number, which may be used without the 

tAs pointed out in Ref. 10 also, there is a 
restriction in the present method on the manner 
in which, for example, At and An go to zero in a 
mesh-refinement procedure. The restriction 
requires that pAt/p(An)^ remain bounded as At and 
An -► p. With this restriction, AT^p(An)^/p^ 
0(An)^* This limitation on Ax is a nuisance 
which one would like to avoid. However, this is 
the price paid for using a simple bidiagonal 
Inversion in place of the more complex 
trldiagonal-inversion procedures for the viscous 
terms . 


CN 


. n 7s 1 

< 0.75 L sJ 




where Cg has also been used in Eq. (21) and can 
have a value of unity or less. The local minimum 
time step near the surface (Atj^)jj„Q, used in 
the above relation, depends on the mesh size 
employed there. Therefore, in this case where 
the shock is treated explicitly and the mesh size 
near the surface is established by requiring that 
the mesh Reynolds number be unity, there would be 
an upper limit on the value of CN which may be 
used with the present method. However, this is 
not a limitation of the method. If the shock 
boundary can be treated implicitly, values larger 
than the limiting value of the Courant number 
indicated here may be used . In fact , the present 
method is unconditionally stable if the flux 

boundary condition 151 6U (see paragraph following 
Eq. (32)) can be evaluated implicitly* It 
represents the implicit part of the boundary 
conditions. Its evaluation by such means as 
lagging in time, etc. will limit the stability of 
the present method to smaller Courant numbers as 
experienced in Refs. 6 and 7. 


Figure 2 gives various time steps which may 
be employed at a given body station. The time 
step shown by curve 2 has been used for CN < 

1 (employing the explicit method^*), whereas "curve 
4 has been employed with the implicit analog^ for 
CN>1. Time step Atj shown by curve 2 is defined 
as 


At^ = (1 - 0.0025 j)A4; j - 1, 2, ..., 100 

where Atj^ is the local minimum time step shown 
by curve 1 and j is the mesh point counter with a 
value of 1 at the wall (n«0) and a value of 100 
just behind the shock. If a very large value of 
the Courant number is used with the implicit 
analogue , a time step shown by curve 5 would 
result and cause the calculations to go implicit 
from the wall to the shock. The global minimum 
curve 3 and the fully-impllcit curve 5 have not 
been used in the present work and are Included 
for illustration only. It becomes clear from 
this figure that, for finer mesh resolution near 
the surface, larger values of the Courant number 
can be used without needing to specify the 
shock-boundary condition implicitly. 

A 101x15 mesh size has been used in the 
present computations with 101 mesh points in the 
direction normal to the surface. The mesh points 
along the body were evenly spaced at A^(=As) 
values of 0.1963 for problem 1 and at A^ values 
of 0.1597 for problem II. The solution Is 
considered as converged to the steady-state value 
when the following criterion is satisfied: 

where e* « 0(10”^) and C|^ is the 

nondimensional heat-transfer coefficient defined 

as 
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The convergence test given here is for CN*10 and 
is made every 50 time steps; for CN«1 it is made 
every 500 time steps, etc. 

During the analysis it was found that the 
quality and stability of the solutions improved 
if the order of differencing along the streamwise 
direction is reversed from one time step to the 
next time step in both predictor and corrector 
steps as given in Eqs. (13) and (14). 

The solution algorithm outlined in this work 
includes the treatment of the source term in the 
implicit parts of the method. However, in the 
results obtained here, the source term in 
implicit parts of Eqs. (13) and (14) was 

neglected, implying that Ax(9^/3U) and 

AtOQb/sS)* etc. were small enough. The 
observation of Ref. 8 in this regard appears to 
be true for the two problems treated here. The 
physics of the problem is contained in the 
explicit parts of Eqs. (13) and (14) and the 
source term is retained there. Dropping the 
source term from the numerics contained in the 
implicit parts of these equations does not seem 
to affect the results appreciably. The inclusion 
of the source term in the implicit part as 
outlined here is likely to increase the 
computational time per time step by an estimated 
5 to 10 percent. 

The computed results for problem I are given 
in Figs. 3 to 5. These figures contain results 
for Courant number (CN) of 1 and 5 with the mesh 
refinement parameter 3 having a value of 1.1. 

The CN=1 results are those of Ref. 12 and have 
been recomputed here. The two factors outlined 
earlier prevented obtaining results for higher 
values of CN. The requirement of keeping the 
mesh Reynolds number around unity did not allow 
further mesh refinement, whereas the specifica- 
tion of the shock boundary condition explicitly 
with the given amount of mesh refinement did not 
permit the use of higher values of CN for this 
problem. 

Figure 3 shows the shock stand-off distance 
along the body surface. There is a very good 
agreement between the two values predicted by 
using CN=1 and 5 over the entire body surface. 

The velocity and temperature profiles of Fig. 4 
show a similar agreement , the comparison being 
somewhat superior over the spherical portion 
(Fig. 4(a)). The comparison between the 
predicted values of skin-friction and heat- 
transfer coefficients and surface pressure for 
CN=1 and 5 is given in Fig. 5. Once again, the 
various distributions compare very well at the 
two CN values. 

Figures 6 through 8 contain results for 
sample problem II. As pointed out earlier, the 
reference Reynolds number in this case is almost 
an order of magnitude larger than the one for 
problem I. Accordingly, a finer mesh (3=1.02) 
can be employed here, still keeping the mesh 
Reynolds number of order unity. This allows the 
use of a higher value of Courant number without 
going implicit all the way up to the shock. 


Thus, results in this case have been obtained for 
the value of Courant number as large as 15. The 
shock stand-off distance of Fig. 6, as well as 
the velocity and temperature profiles of Fig. 7, 
compare quite well at various CN values. Once 
again, the comparison is superior over the 
spherical portion (s»0.32) of the body. Figure 8 
also shows good agreement between the distribu- 
tions of skin-friction and heat-transfer coeffic- 
ients and surface pressure employing three values 
of the Courant number. 

The results presented here have been 
obtained on the Control Data CYBER 203 vector- 
processing computer. The explicit part of the 
method is fully vectorized. The increase in 
computing time due to the addition of implicit 
part is about 70 percent. Typically it takes 
about 3.0x10 ^ s per mesh point per time step for 
the explicit part using the local minimum time 
step.tt With the addition of the implicit part) 
this time increases to about 5x10”^ s. Two 
factors which influence the computing time in the 
present case are the partial vectorization of the 
implicit part and the extent to which the 
implicit part is called if the shock boundary is 
treated explicitly. However, even with the 
increased computing time per time step, the total 
computing time is reduced significantly. In case 
of problem II, for example, the total computing 
time is about 6 times less at CN=15 as compared 
to the computing time at CN=1 with the mesh 
refinement parameter 3 having a value of 1.02. 

For a higher mesh resolution near the surface, a 
larger value of the Courant number may be used 
while still treating the shock explicitly. This 
would reduce the total computing time even 
further when compared to the computing time at 
CN = 1. The mesh refinement, however, is done 
only to the point where the mesh Reynolds number 
would be of order unity or so for the reasons 
explained earlier. 


Concluding Remarks 

In the analysis presented here, the 
recently reported implicit analog of MacCormack*s 
earlier widely used explicit method has been 
extended to external axisymmetric laminar flows 
with strong entropy gradients. The details of 
the ‘‘numerics" of the implicit part are obtained 
in a body-oriented coordinate system with a 
moving outer (shock) boundary during the 
transient part of the solutions. The implicit 
analog is unconditionally stable if the boundary 
conditions are also specified fully Implicitly 
(i.e. without the lag of 1/2 or 1 time step). In 
this work, the inner (wall) boundary for the 
implicit part was treated through the reflection 
concept as suggested in the original presentation 
of the implicit analog by one of the present 
authors (RWM) • The outer (shock) boundary was 
treated explicitly in order to avoid the 
specification of this boundary condition 
Implicitly and to keep the computational 
algorithm simple. Thus, the present results do 
not contain any approximation about the treatment 
of the boundary condition which may affect the 

TTThe same computing time per mesh point per time 
step for the explicit part is required if the 
global minimum time step is employed. However, 
the total computing time in this case is 
increased by a factor of about 2.5. 


9 



stability of the implicit analog other than the 
Courant number limitation. The limiting values 
of the Courant number are provided when the shock 
boundary is treated explicitly. The solution 
algorithm outlined includes the treatment of the 
source term present in a weakly conservative 
system of equations. 

For the results presented here, the Courant 
number in the coordinate direction along the body 
surface was 0(1), but the Courant number in the 
direction normal to the body, because of the fine 
mesh-point spacing needed to resolve the viscous 
boundary layer at the body surface, varied from 5 
to 15 depending on flow Reynolds number. A 
detailed comparison of various results obtained 
with different values of the Courant number (CN) 
shows that the accuracy of predictions is quite 
good at higher values of CN. There is a signifi- 
cant saving in the overall computing time even 


though the computing time per time step increases 
by about 70 percent with the inclusion of the 
implicit part. The code developed for the 
implicit analogue on the Control Data CYBER 203 
computer was essentially an "addition*' to the 
earlier explicit code. The coding of the 
Implicit part is not affected much with the 
increase in complexity of the explicit part. 

This may have special appeal in the analysis of 
3-D flow problems. 
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and the transformation Jacobians J and J are 
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Appendix B; Jacobian Matrices 1 and 5 
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The Jacobian matrices A and S are obtained from 


«t 


■ 

0 

n Qt * 8 * 

5 + (i-e*)ic3 

+U5 

ilf 
X ^8 

-vu + a*3*Cj^ + 

X«8 - 

u + (1-S')v5jj 


(-SHt)l(l-e’)a’ +p-] 

c2 


(e’+i)S 


+ u3 - u3 *u 

+ - vg'u 

-3*5t 



Is 

X 

% 

0 

“.0*3* 
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V + (1-3*) j tig 

3 *v . 

- — ris + un„ 
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CD 

“VV + a*3*Hjj + vn^ 

X'<s ■ 

V + (1-3*)vti^ 


(-v+nt)((l-3')a' +1^] 

/C^ . IX 

<3^ + “ > T 


(3'+l)v 


+ u3 *n^ “ u3 *v 

+ v3n^ - v3*v 

-3X 


(Bl) 


(B2) 


where c = VYP/p is the speed of sound, a* = (u^+v^)/2, and 3* » (Y“1)« 


Appendix C; Matrices Sp , , and 

The matrices , and S^^are given by 


1 _ “'6' 
c2 

u3* 

c2 

v3 ' 
c2 

r 

’Ji^o 

1 

c2 « 

c('3* -|- (u-C^) 

4^3 -«e* 

c2 

3* 

(ucC^- vc5g/X) 

pd 

^ E 

pd ^n 

Xpd ^s 

0 

a'3' +|- (v-n^) 

-4«8 

, c2 

-r-«n 

3* 


(Cl) 
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1 

2c^ 


0 








^ «n> 


C-UC_ + 7- CJ 


(vCTig/X - UCIIjj) 


c2 « 

a'6' (v-n^) 




c 

Xpg ''s 


-jr,„-ve* 





r C^ . 

-xi%) 

2c^ 

CM 

O 

CM 

1 

(V +4'ln) 

(V - T,^) 

CM 

2c^ 


11, -i„ 

* * “ » 2c" 2c" 


If it is assumed that the outer boundary (which is shock in the present case 
and moves with time during the transient part of the solutions) is fixed, 
for example for the internal flow problems, then nt are 

identically zero and Eqs, (Cl) through (C4) reduce to those obtained in 
Ref. 7 for a flat surface (A“l)* It may be mentioned that for obtaining the 
expressions provided here certain columns and rows of Eqs. (12) and (13) of 
Ref. 7 need to be multiplied and divided by (pay2) and a change in their 
order is required. 
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Fig. 2 Definitions of various time steps at the 
stagnation point for Problem II and 
Courant number (CN) of 5 and mesh 
refinement parameter 0 » 1.02. 
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Fig. 3 Shock stand-off distance for two values 
of the Gourant number with B * 1.1, 
Problem I. 
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Fig. 4 Concluded. 
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(a) SKIN-FRICTION COEFFICIENT 


Fig. 5 Comparison of wall quantities with two 

values of the Courant number for B * 1^1, 
Problem I. 


Fig. 4 Comparison of velocity and temperature 

profile for two values of Courant number 
over the spherical portion of the body 
for B “ 1*1» Problem I. 
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Fig • 5 Continued • 






Fig, 6 Shock stand-off distance obtained with 
different values of the Courant number 
for 3 = 1.02, Problem II. 
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(a) TANGENTIAL VELOCITY 

7 Comparison of flowflald profiles over the 
spherical (s - 0.32) and conical flank 
(s « 1.76) portions of the body for 
6 - 1.02, Problem II. 
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Fig. 8 Comparison of wall quantities with 

different values of the Courant number 
(CN) for 3 - 1.02, Problem II. 
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